# Run on R 3.3.0 GUI 1.68 Mavericks build (7202)
# June 3, 2016

## This code runs, in this order: 
 # 1. HLM model in column 6 of Table 1. 
 # 2. Figure 3.
 # 3. HLM model in column 2 of Table 2. 
 # 4. Figure 4. 


#############################################
# 2-LEVEL HIEARCHICAL MODELS 
############################################



library(lme4)
library(SASmixed)
library(foreign)
library(texreg)
library(arm)
library(lmerTest)

#### Set Path
#setwd(...)

dataindustry <- read.dta("industry_level_data_for_HLM_models.dta")

# Retrieve individual predictors
country <- dataindustry$country
idc <- dataindustry$idc
wtaxcomp <- dataindustry$wtaxcomp
tariff <- dataindustry$tariff
region <- dataindustry$region
avgcompetitors <- dataindustry$avgcompetitors
lnavgexp <- dataindustry$lnavgexport
obsolete <- dataindustry$obsolete
lntariff <- dataindustry$lntariff
lnavgexp <- dataindustry$lnavgexport
lnavgage <- dataindustry$lnavgage
lntotlabor <- dataindustry$lntotlabor
gov_owner_share <-dataindustry$gov_owner_share
avgexp <- dataindustry$avgexport
ln_tot_pop <- dataindustry$ln_tot_pop
mining  <- dataindustry$mining
free_media <- dataindustry$free_media

n <- length(tariff)

# Retrieve country-level predictors
datacountry <- read.dta("country_level_data_for_HLM_models.dta")
lowfc <- datacountry$lowfc
idc <- datacountry$idc
		
# Create country index variable 
country.name <- as.vector(dataindustry$country)
uniq <- unique(country.name)
J <- length(uniq)
country <- rep (NA, J)
for (i in 1:J){
  	country[country.name==uniq[i]] <- i
	}

# Assign Taxstaff value to Country i
lowfc.full <- lowfc[country]
idc.full <- idc[country]

# Interaction
inter <- tariff*obsolete
	

########################################
# HLM Model, Column 6, TABLE 1 
########################################

C6 <- lmer(formula = wtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 + tariff + obsolete + inter | country), verbose = TRUE)

## P-values calculated based on Satterthwate's approximations [lmerTest]
summary(C6)

## Anovat Test: Are random effects statistically significant?

C6.null <- lmer(formula = wtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 | country), verbose = TRUE)

C6.full <- lmer(formula = wtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 + tariff + obsolete + inter | country), verbose = TRUE)

anova(C6.null, C6.full)


######################################################
# FIGURE 3. Estimates drawn from Column 6 in TABLE 1
######################################################

# Obsolete X Tariff
b3.hat.C6 <- fixef(C6)[13] + fixef(C6)[17]*lowfc + ranef(C6)$country[,4]
b3.se.C6 <- se.ranef(C6)$country[,4]

lower.b3 <- b3.hat.C6 - b3.se.C6
upper.b3 <- b3.hat.C6 + b3.se.C6

# pdf(".../Figure_3.pdf", width=10,height=7)	
plot (lowfc, b3.hat.C6, cex.lab=1.2, cex.axis=1.1, ylim=range(lower.b3,upper.b3),
xlab="Fiscal Capacity", ylab="Three-way interaction coefficient", pch=20)
# Zero-Line
lines(lowfc, rep(0, length(lowfc)), type = "l", col = "grey90")
# Grid for background
grid(col = "gray80")
mtext("High",side=1, line=3, at=-2, col="black")
mtext("Low",side=1, line=3, at=0, col="black")
curve (fixef(C6)[13] + fixef(C6)[17]*x, lwd=1, col="black", add=TRUE)
segments (lowfc, lower.b3, lowfc, upper.b3, lwd=.5, col="gray10")		
#dev.off()



########################################
# HLM Model, Column 2, TABLE 2
########################################

# Unweighted Tax Compliance
uwtaxcomp <- dataindustry$uwtaxcomp

C2 <- lmer(formula = uwtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 + tariff + obsolete + inter | country), verbose = TRUE)
## P-values calculated based on Satterthwate's approximations [lmerTest]
summary(C2)


#### Anovat Test: Are random effects statistically significant?

# Level 3 random effects are jointly statistically 
C2.null <- lmer(formula = uwtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 | country), verbose = TRUE)

C2.full <- lmer(formula = uwtaxcomp ~ region + mining + lnavgexp + avgcompetitors + lnavgage + lntotlabor + gov_owner_share + free_media + ln_tot_pop + obsolete + tariff + inter + lowfc.full + tariff:lowfc.full + obsolete:lowfc.full + inter:lowfc.full + (1 + tariff + obsolete + inter | country), verbose = TRUE)

anova(C2.null, C2.full)


######################################################
# FIGURE 4. Estimates drawn from Column 2 in TABLE 2
######################################################

# Obsolete X Tariff
b3.hat.C2 <- fixef(C2)[13] + fixef(C2)[17]*lowfc + ranef(C2)$country[,4]
b3.se.C2 <- se.ranef(C2)$country[,4]

lower.b3 <- b3.hat.C2 - b3.se.C2
upper.b3 <- b3.hat.C2 + b3.se.C2

# pdf(".../Figure_4.pdf", width=10,height=7)	
plot (lowfc, b3.hat.C2, cex.lab=1.2, cex.axis=1.1, ylim=range(lower.b3,upper.b3),
xlab="Fiscal Capacity", ylab="Three-way interaction coefficient", pch=20)
# Zero-Line
lines(lowfc, rep(0, length(lowfc)), type = "l", col = "grey90")
# Grid for background
grid(col = "gray80")
mtext("High",side=1, line=3, at=-2, col="black")
mtext("Low",side=1, line=3, at=0, col="black")
curve (fixef(C2)[13] + fixef(C2)[17]*x, lwd=1, col="black", add=TRUE)
segments (lowfc, lower.b3, lowfc, upper.b3, lwd=.5, col="gray10")		
#dev.off()




